Optimal Dynamical Range of Excitable Networks at Criticality 
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, A recurrent idea in the study of complex systems is that optimal information processing is to 

be found near bifurcation points or phase transitions. However, this heuristic hypothesis has few 
(if any) concrete realizations where a standard and biologically relevant quantity is optimized at 
^-^1 criticality. Here we give a clear example of such a phenomenon: a network of excitable elements 

. has its sensitivity and dynamic range maximized at the critical point of a non-equilibrium phase 

transition. Our results are compatible with the essential role of gap junctions in olfactory glomeruli 
and retinal ganglionar cell output. Synchronization and global oscillations also appear in the net- 
1/^ . work dynamics. We propose that the main functional role of electrical coupling is to provide an 

enhancement of dynamic range, therefore allowing the coding of information spanning several orders 
of magnitude. The mechanism could provide a microscopic neural basis for psychophysical laws. 
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' Psychophysics is probably the first experimental area in neuropsychology, having been founded by physiologists 
QT'^ (Helmholtz and Weber) and physicists (Fechner, Plateau, Maxwell and Mach) in the middle of the nineteenth cen- 
tury Its aim is to study how physical stimuli transduce into psychological sensation, probably the most basic 
' mind-brain problem. How to relate such high-level psychological phenomena to low-level neurophysiology is a task 
^ ' for the twenty-first century brain sciences. 

I As the intensity of physical stimuli (light, sound, pressure, odorant concentration etc) varies by several orders of 
' magnitude, psychophysical laws must have a large dynamical range. Although this has been extensively verified 
' experimentally at the psychological ^1;] and neural i4|| level, little work has been done regarding the mechanism 
that produces such psychophysical laws However, two main ideas seem to be consensual: (1) nonlinear 

\^ transduction must be done at the sensory periphery to prevent early saturation and (2) the broad dynamic range 
is often a collective phenomenon [l, i5|, i6|] , because single cells usually respond in a linear saturating way with small 
"q ■ ranges [9l.[loj. 

• ^ ' Two nonlinear transfer functions have been widely used to fit experimental data, both in psychophysics and in neural 
1^ response: the logarithm function F{S) — Clog 5 (Weber-Fechner law) and the power law function F{S) ~ CS"^ 
O^" (Stevens Law), where S is the stimulus level, C is a constant and m is the Stevens exponent. Later, other functions 
^ ' have been proposed to fit data with more extended input range, and to account for sensory saturation, in particular 
the Hill function F{S) = F^axS"^ / {S"^ -I- S*™) where F^ax is the saturation response and Sq is the input level for 
half-maximum response. Notice that, because both Hill and Stevens functions have a power law regime, it is natural 



• to denote the exponents by the same parameter m. 

Some authors have tried to derive such phenomenological laws from the structure of natural signals [11]. This 
type of work may furnish an evolutionary motivation for biological organisms to implement or approximate such laws. 
However, there is no consensual view on this theoretical aspect, and it does not provide a neural basis for implementing 
the psychophysical laws. 

In contrast, our statistical physics approach to neural psychophysics 0, d, shows how Hill-like transfer 
functions may arise in biological excitable systems: they are not put into the models by hand, nor mathematically 
derived from a priori reasoning, but appear as a cooperative effect in a network of excitable elements. At diflterente 
biological levels, these elements may be interpreted as whole neurons, excitable dendrites, axons, or other subcellular 
excitable units. We show that a network of excitable elements, each with small dynamic range, presents a collective 
response with broad dynamic range and high sensitivity. Even more interesting, we find that the dynamic range 
is maximized if the spontaneous activity of the network corresponds to a critical process. This is compatible with 
recent findings of critical avalanches in in vitro neural networks 13, l3| and provides a clearcut example of optimal 
information processing at criticality 15, 3, l3|- The model also has other dynamical features and permits us to make 
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a testable prediction. 

In previous work we have introduced the idea that excitable waves in active media provide a mechanism for strong 
nonlinear amplification with large dynamic range. This has been shown in simulations with one-dimensional [y| 
and two-dimensional p?] deterministic cellular automaton models, as well as with one-dimensional networks of coupled 
maps and Hodgkin-Huxley elements Q . Analytical results have recently been obtained for the one-dimensional cellular 
automaton model under the two-site mean-field approximation [l^ . In this paper, we study a very different system 
where the activity propagation is stochastic and the electrical synapses form a random network. This latter case 
seems to be a more realistic topology for, say, the olfactory intraglomerular network of excitable dendrites coupled 
by gap junctions 3, 3 2^- We find a whole new phenomenology related to a non-equilibrium phase transition to 
re-entrant activity present in this model. 

In the present model, each excitable element i — 1, . . . , N has n states: Si = is the resting state, Si = 1 corresponds 
to excitation and the remaining = 2, . . . , n — 1 are refractory states. There are two ways for the ith element to 
go from state Si = to 1: (1) owing to an external stimulus, modelled here by a Poisson process with rate r (which 
implies a transition with probability A = 1 — exp(— rAt) per time step); (2) with probability , owing to a neighbour 
j being in the excited state in the previous time step. Time is discrete (we assume At — 1 ms) and the dynamics, 
after excitation, is deterministic: if Sj = 1, then in the next time step its state changes to Si = 2 and so on until the 
state Si = n — 1 leads to the Si — resting state, so the element is a cyclic cellular automaton [2l[. The Poisson 
rate r will be assumed to be proportional to the stimulus level S (for example, the odorant concentration in olfactory 
processing). Notice that each element receives external signals independently, that is, we have a Poisson process for 
each element (modelling the arrival of axonal inputs from different receptor neurons). 

The network with N elements is an Erdos-Renyi undirected random graph, with NK/2 links being assigned to 
randomly chosen pairs of elements. This produces an average connectivity K where each element i (i = 1, . . . , N) is 
randomly connected to Ki neighbours. The distribution P{Ki) of neighbours is a binomial distribution with average 
K (see Fig fA). The probability that an active neighbour j excites element i is given by pij, a random variable 
with uniform distribution in the interval [Q,Pmax]- The weights are symmetrical {pij = pji) and are kept fixed 
throughout each simulation ("quenched disorder"). This kind of coupling models electric gap junctions instead of 
chemical synapses because it is fast and bi-directional. However, symmetry is not a necessary ingredient, because 
similar results are obtained in asymmetrical networks. Note that we are not assuming that gap junctions have a 
stochastic dynamics, but only that, because other internal factors and noise are also present, a probabilistic account 
of the excitation process may be more realistic. 

The local branching ratio ctj = X]f ^ Pij corresponds to the average number of excitations created in the next time 
step by the j-th element fl^]. The distribution P{(Jj) of local branching ratios (a bell-shaped distribution with average 
a = {(Ji)) is shown in Fig. lA. The average branching ratio a is the relevant control parameter. In the simulations, 
we set a by choosing Pmax = 2a/K and keeping a < K/2. 

The network instantaneous activity is the density pt of active (s — 1) sites at a given time t. We also define the 
average activity F ~ T^^ X)t=i Pt where T is a large time window (of the order of 10^ time steps). Typical pt curves in 
the absence of stimulus (r = 0) are shown in Fig. IB. Notwithstanding the large variance of P{aj)^ only supercritical 
networks (that is, with ct > CTc = 1) have self-sustained activity, F > 0. As expected, critical networks have a larger 
variance in the distribution of extinction times and present a power law behaviour in the distribution of avalanche 
sizes with exponent 3/2 (not shown), in agreement with findings in biological networks [isj . 

Two types of oscillations are observed in this system. Under sufficiently strong stimulation, all networks present 
transient collective oscillations, with frequencies of the order of the inverse refractory period (Fig. IC). They are a 
simple consequence of the excitable dynamics and the sudden synchronous activation by stimulus initiation. This 
transient behaviour is reminiscent of oscillations widely observed in experiments [2^ 1 . Networks with a > aosc > o'c 
also present self-sustained oscillations in the absence of stimulus (Fig. ID), where aosc is a bifurcation threshold. The 
frequency depends on the network parameters, but remain in the gamma range 30 — 60 Hz. The oscillations are similar 
to reentrant activity found in spatially extended models of electrically coupled networks (23| , from which analytical 
techniques for calculating the frequency could perhaps be borrowed [24| . We remark that the large time window used 
in the definition of F is only chosen for convenience, as it can be seen from Figs. 1C,D and 200 ms is sufficient to give 
reliable averages. 

As a function of the stimulus intensity r, networks have a minimum response Fq (= for the subcritical and critical 
cases) and a maximum response Fmax- We define the dynamic range A = 10 log(ro.9/7'o.i) as the stimulus interval 
(measured in dB) where variations in r can be robustly coded by variations in F ^ discarding stimuli which are too 
weak to be distinguished from i^o or too close to saturation. The range [?'o.i,ro.9] is found from its corresponding 
response interval [Fo.i,Fo.g], where F^ ^ Fq + x{Fmax — Fq) (see Fig. 2C). This choice of a 10%-90% interval is 
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FIG. 1: Network characterization and density of active sites. Simulations with N — 10^ sites, K — 10 and n = 10 states. 
(A) Probabihty density function of local branching ratio and (inset) connectivity ("degree") distribution; (B) Instantaneous 
density of active sites for subcritical (black) , critical (red) and supercritical (blue) branching parameters as functions of time 
(three different runs for each case); (C), (D) Instantaneous density (of all sites; upper panels) and raster plot (of 10^ randomly 
chosen sites; lower panels) in response to a square pulse of stimulus (r = 0.5 ms~^ for 100 ms < t < 300 ms, null otherwise) 
for critical (C) and supercritical (D) branching parameters. 



arbitrary, but is standard in the literature and does not affect our results. 

As can be seen in Figs. 2A and 2B, the response curves F{r) of the networks present a strong enhancement of 
dynamic range compared with the uncoupled case tr = 0. In the subcritical regime, sensitivity is enlarged because 
weak stimuli are amplified due to activity propagation among neighbours. As a result, the dynamic range A(cr) 
increases monotonically with a. In the supercritical regime, the spontaneous activity Fq masks the presence of weak 
stimuli, therefore A(cr) decreases. The optimal regime occurs precisely at the critical point (see Fig. 2D). This is a 
new and important result, because it is perhaps the first clear example of signal processing optimization at a phase 
transition, making use of a standard and easily measurable performance index. 

The curves F{r) could be fitted by a Hill function, but are not exactly Hill. The theoretical curves in Fig. 2 are 
obtained from a simple mean-field calculation (see below) that provides a very good fit to the simulation data without 
free parameters, and correctly predicts the exponents governing the low-stimulus response F cx r"*. An important 
point is that the Stevens-Hill exponent m changes from m = 1 in the subcritical regime to m — 0.5 at criticality. If we 
assume that biological networks work in the optimal regime, the critical value m — 0.5 suggests how exponents less 
than one could emerge in psychophysics [if and neural responses f2, 4]. We note that apparent exponents between 
0.5 and 1.0 are observed [l2] if finite size effects are present, that is, if N is small. 

In the simple mean- field approximation, we take Ki — K and pij as the average value a/K. The probability pt 
that an inactive site at time t will be activated in the next time step by at least one of its K neighbours (a fraction 
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FIG. 2: Response curves and dynamic range. Points represent simulation results with A*" = IC sites, K — 10, n = 5 states and 
T = 10^ ms, whereas lines correspond to the mean-field model described in the text. (A) Response curves (mean firing rate vs. 
stimulus rate) from a = to cr = 2 (in intervals of 0.2). Line segments are power laws F oc with m — 1 (subcritical) and 
m = 1/2 (critical). Inset: spontaneous activity -Fb vs. branching ratio a. (B) The same as in (A), but with a linear vertical 
scale. (C) Response curve for a = 1.2 and relevant parameters for calculating the dynamic range A. (D) Dynamic range vs. 
branching ratio is optimized at the critical point a — 1. 



Ft of which is currently active) is simply = 1 — (1 — crFt/K)^ . This leads to the following mean-field map: 

Ft+i ^ PtiO)X + Ptmi - X)pt , (1) 

where 



Pt{0) ^ 1 ~ {n ~ l)Ft (2) 

is the approximate probability of finding a site in the resting state and A(r) = 1 — exp(— rAt). The first term on the 
right-hand side of equation [T] corresponds to activation due to external input, and the second term corresponds to 
activation due to neighbour propagation. In the stationary regime, the response function F(r) is given by the solution 
of 

F = (1 - (n - 1)^^) [1 - (1 - aF/K f (1 - A(r))] . (3) 

Note that equation [2] is only exact in the stationary state, but the resulting map is consistent because the result of 
its iterations coincides with the solution of equation [3] (in particular, F < 1/n, so that P(0) > 0). 



As is usual in the statistical physics of phase transitions [21|, we analyse two limits: the critical behavior without 



an external field and the effect of a vanishing field at the critical point. In the absence of external stimulus (A = 0), in 
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the limit cr — > Cc = 1 (-P ^ 0), the order parameter behaviour is F{a) ~ (cr — 1)/^ where C = (n — l) + {K — 1)/2K , 
which gives, from the definition F(a) oc (cr — (Tc)'^, a critical exponent = 1. On the other hand, at the critical point 
cr = CT = 1 we have F{r) « \/r/C which gives the Stevens exponent m = 1/2. Notice that the rate r plays the role 



of an external field h and, from the usual statistical-physics definition 21| F oc h}/^*^ we get the critical exponent 
5h = 2. As expected, these are the classical mean-field exponents for branching processes, which seem to be valid 
for our model even with the presence of quenched disorder. An important conceptual point is that the Stevens-Hill 
exponent is indeed found to be indeed a statistical-physics critical exponent m — l/Sh- 

Optimization of dynamic range at criticality is very robust; the shape of the plot in Fig. 2D does not depend on 
parameters such as the average number of neighbours K or refractory period n. The random-network case seems 
to be a lower bound for the dynamic-range enhancement. In future works, we will report that A(cr) is even more 
enhanced in low-dimensional networks, and presents non-rational Stevens exponents. Perhaps a better compromise 
between larger dynamic range and biological realism would be a small-world or scale-free network. 

As an example of the robustness of the results, if we put Hodgkin-Huxley elements with coupling with values = 
or 1, which would correspond to a deterministic case with the presence of absence of bonds, we obtain a A{p) curve, 
where p is the probability that a bond exists between the elements. This A{p) curve is very similar to the A(cr) 
curve of the present model, that is, A(p) has a strong peak at the percolation phase transition p — pc (to be reported 
elsewhere). However, the present model has the virtue of enabling analytical results that provide a benchmark for 
the performance of networks with other topologies. 

Now we discuss the possible relevance of our results to biological sensory processing. Recent findings show that 
projection cells in sensory systems are coupled via dendro-dendritic electrical synapses, for example a-ganglionar cells 
in the retina [i^ [i^ and mitral cells in the olfactory bulb 13, 3, 23] ■ In most of these findings the electrical coupling 
is mediated by connexins, but pannexins could also be present, and could even be more important than connexins 
for providing electrical coupling between excitatory cells 27|. However, the functional role of this electrical coupling 
is largely unknown. We propose that the electrically coupled dendritic trees in these systems form an excitable 
network (where each element of our model represents an excitable dendritic patch). Our results are consistent with 
the reduction in sensitivity, dynamical range and synchronization recently observed in retinal ganglion cell response 
of connexin-36 knockout mice . 

In the case of olfactory system, we identify the excitable random network with the dendro-dendritic network in 
the glomeruli ,18], and each element is interpreted as an active dendritic compartment containing ion channels. It is 
known that relevant electrical coupling between mitral cells is done at the glomerular level 13, 20 1 because only cells 
that have their apical dendrite tufts in the same glomerulus show synchronized activity. In connexin-36 knockout 
mice, the synchronized activity of mitral cells is absent. 

Our hypothesis could be tested in the following way. The dynamic range of glomeruli is of the order of 30 dB as 
measured recently ^,'3] (in contrast to 10 dB of single olfactory receptor neurons). We predict that in connexin-36 
knockout mice, this dynamic range shall be strongly reduced. Of course, for a decisive test, we need to examine if 
other electrical synapses based on connexin- 45 [H and pannexins ^ are irrelevant in olfactory glomeruli. 

The text-book account of large dynamic range in intensity coding is the recruitment model or its variants [5] where 
different elements, with diverse activation thresholds, are sequentially recruited. Our mechanism is not incompatible 
with this scenario, and we expect that threshold variability in our model enlarges the dynamic range even more. 
However, we must remember that recruitment is a linear mechanism where, to explain each order of magnitude in 
dynamic range, we must postulate a corresponding order of magnitude in activation thresholds, which is not a plausible 
assumption [51]. 

We may ask how the network could self-organize to the critical point (Tc = 1. It is not hard to conceive that 
homeostatic mechanisms, acting on the number and conductance of gap junctions, could tune the system. It is well 
known that extensive pruning of gap junctions occurs during development and maturation. This could represent the 
initial self-organization process towards criticality. Spontaneous activity in the absence of input furnishes a signature 
of supercriticality that could be used as a feedback signal to control the system. 

This self-tuning criticality has recently been proposed in a model of sound nonlinear amplification by Hopf oscillators 
in the cochlea [30[. In the cochlear model, it is assumed that each oscillator is poised at its Hopf bifurcation point, 
producing enhanced sensitivity and enlarged dynamic range. The principal difference from our mechanism is that it 
is a model for individual cells, that is, it is not based on a collective phenomenon. This implies that the dynamic- 
range exponent must be classical (rational). In our case, the rational exponent m — 1/2 is a particularity of the 
random network, and, similarly to other statistical-mechanics models, more-structured network topologies may have 
non-classical exponents. 

Although in this work we restricted our attention to sensory processing, we note that the dynamic range of more 
central networks could also be improved by the same mechanism. Similarly to excitatory networks, inhibitory networks 
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must also work robustly in the presence of large variations in input. So, the presence of electrical synapses in cortical 
inhibitory networks [31l] could reflect the same principle. 

The mechanism for amplified nonlinear response due to wave creation and annihilation is a basic property of 
excitable media. We found in this work that, if active media are tuned at the critical point of activity propagation, 
the response is optimized. We proposed that this principle is present in electrically coupled excitable dendritic networks 
in projection neurons of sensory systems and is a generative mechanism for psychophysical laws. This computational 
principle based in critical activity could also be present in other brain regions and could be implemented in artificial 
sensors by using excitable media as detectors. 

This research is supported by CNPq, FACEPE, CAPES and PRONEX. The authors are grateful for discussions 
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